Подгружаем датасет
led <- read_rds("./life_expectancy_data.RDS")
Смотрим на переменные
str(led)
## Classes 'data.table' and 'data.frame': 195 obs. of 23 variables:
## $ Country : chr "Afghanistan" "Albania" "Algeria" "Angola" ...
## $ Year : int 2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
## $ Gender : chr "Female" "Female" "Female" "Female" ...
## $ Life expectancy : num 66.4 80.2 78.1 64 78.1 ...
## $ Unemployment : num 14.06 11.32 18.63 7.84 8.26 ...
## $ Infant Mortality : num 42.9 7.7 18.6 44.5 5.1 ...
## $ GDP : num 18799450743 15400242875 171767403748 89417190341 1687533333 ...
## $ GNI : num 19098305631 15198661275 167568133680 81900890781 1581477852 ...
## $ Clean fuels and cooking technologies : num 36 80.7 99.3 49.6 100 ...
## $ Per Capita : num 494 5396 3990 2810 17377 ...
## $ Mortality caused by road traffic injury: num 15.9 11.7 20.9 26.1 0 ...
## $ Tuberculosis Incidence : num 189 16 61 351 0 29 26 2.2 6.9 6 ...
## $ DPT Immunization : num 66 99 91 57 95 ...
## $ HepB3 Immunization : num 66 99 91 53 99 ...
## $ Measles Immunization : num 64 95 80 51 93 ...
## $ Hospital beds : num 0.432 3.052 1.8 0.8 2.581 ...
## $ Basic sanitation services : num 49 99.2 86.1 51.4 85.5 ...
## $ Tuberculosis treatment : num 91 88 86 69 72.3 ...
## $ Urban population : num 25.8 61.2 73.2 66.2 24.5 ...
## $ Rural population : num 74.2 38.8 26.8 33.8 75.5 ...
## $ Non-communicable Mortality : num 36.2 6 12.8 19.4 17.6 ...
## $ Sucide Rate : num 3.6 2.7 1.8 2.3 0.8 ...
## $ continent : Factor w/ 5 levels "Africa","Americas",..: 3 4 1 1 2 2 4 2 5 4 ...
## - attr(*, ".internal.selfref")=<externalptr>
## - attr(*, "sorted")= chr "Country"
summary(led)
## Country Year Gender Life expectancy
## Length:195 Min. :2019 Length:195 Min. :55.49
## Class :character 1st Qu.:2019 Class :character 1st Qu.:70.02
## Mode :character Median :2019 Mode :character Median :77.55
## Mean :2019 Mean :75.52
## 3rd Qu.:2019 3rd Qu.:80.95
## Max. :2019 Max. :88.10
## Unemployment Infant Mortality GDP
## Min. : 0.178 Min. : 1.40 Min. : 188391771
## 1st Qu.: 3.735 1st Qu.: 5.35 1st Qu.: 11170154590
## Median : 5.960 Median :13.50 Median : 39670977333
## Mean : 8.597 Mean :19.61 Mean : 465995962784
## 3rd Qu.:10.958 3rd Qu.:30.23 3rd Qu.: 247555140295
## Max. :36.442 Max. :75.80 Max. :21433224697000
## GNI Clean fuels and cooking technologies
## Min. : 375392118 Min. : 0.00
## 1st Qu.: 10943629972 1st Qu.: 34.50
## Median : 40087702395 Median : 80.70
## Mean : 486438131433 Mean : 65.98
## 3rd Qu.: 245717112337 3rd Qu.:100.00
## Max. :21708650000000 Max. :100.00
## Per Capita Mortality caused by road traffic injury
## Min. : 228.2 Min. : 0.00
## 1st Qu.: 2165.3 1st Qu.: 8.20
## Median : 6624.8 Median :16.00
## Mean : 16821.0 Mean :17.06
## 3rd Qu.: 19439.7 3rd Qu.:24.00
## Max. :175813.9 Max. :64.60
## Tuberculosis Incidence DPT Immunization HepB3 Immunization
## Min. : 0.0 Min. :35.00 Min. :35.00
## 1st Qu.: 12.0 1st Qu.:85.69 1st Qu.:81.31
## Median : 46.0 Median :92.00 Median :91.00
## Mean :103.8 Mean :87.99 Mean :86.76
## 3rd Qu.:138.5 3rd Qu.:97.00 3rd Qu.:96.00
## Max. :654.0 Max. :99.00 Max. :99.00
## Measles Immunization Hospital beds Basic sanitation services
## Min. :37.00 Min. : 0.200 Min. : 8.632
## 1st Qu.:84.85 1st Qu.: 1.301 1st Qu.: 62.919
## Median :92.00 Median : 2.570 Median : 91.144
## Mean :87.31 Mean : 2.997 Mean : 77.380
## 3rd Qu.:96.50 3rd Qu.: 3.773 3rd Qu.: 98.582
## Max. :99.00 Max. :13.710 Max. :100.000
## Tuberculosis treatment Urban population Rural population
## Min. : 0.00 Min. : 13.25 Min. : 0.00
## 1st Qu.: 73.00 1st Qu.: 41.92 1st Qu.:21.98
## Median : 82.00 Median : 58.76 Median :41.24
## Mean : 77.57 Mean : 59.12 Mean :40.88
## 3rd Qu.: 88.00 3rd Qu.: 78.02 3rd Qu.:58.08
## Max. :100.00 Max. :100.00 Max. :86.75
## Non-communicable Mortality Sucide Rate continent
## Min. : 4.40 Min. : 0.300 Africa :52
## 1st Qu.:11.85 1st Qu.: 2.050 Americas:38
## Median :17.20 Median : 3.500 Asia :42
## Mean :17.05 Mean : 4.802 Europe :48
## 3rd Qu.:22.10 3rd Qu.: 6.600 Oceania :15
## Max. :43.70 Max. :30.100
Все выглядит пристойно.
Оценим пропущенные значения.
led %>%
is.na() %>%
sum()
## [1] 0
Их нет, и это прекрасно.
led %>%
select(where(is.numeric))
## Year Life expectancy Unemployment Infant Mortality GDP
## 1: 2019 66.388 14.065000 42.90000 18799450743
## 2: 2019 80.201 11.322000 7.70000 15400242875
## 3: 2019 78.133 18.629999 18.60000 171767403748
## 4: 2019 64.039 7.835000 44.50000 89417190341
## 5: 2019 78.104 8.256285 5.10000 1687533333
## ---
## 191: 2019 79.523 2.008000 14.50000 261921244843
## 192: 2019 83.100 13.543000 30.23477 4068000000
## 193: 2019 67.826 25.027000 42.20000 24749812134
## 194: 2019 66.891 13.089000 38.90000 23308667781
## 195: 2019 62.899 5.248000 33.80000 19284289739
## GNI Clean fuels and cooking technologies Per Capita
## 1: 19098305631 36.00000 494.1793
## 2: 15198661275 80.70000 5395.6595
## 3: 167568133680 99.30000 3989.6683
## 4: 81900890781 49.60000 2809.6261
## 5: 1581477852 100.00000 17376.6497
## ---
## 191: 245123373262 64.70000 2715.2760
## 192: 2042053577026 60.59302 38136.6658
## 193: 23509163481 60.90000 1055.0502
## 194: 22904080807 15.70000 1305.0010
## 195: 18896621255 30.00000 1316.7407
## Mortality caused by road traffic injury Tuberculosis Incidence
## 1: 15.90000 189.0000
## 2: 11.70000 16.0000
## 3: 20.90000 61.0000
## 4: 26.10000 351.0000
## 5: 0.00000 0.0000
## ---
## 191: 30.60000 176.0000
## 192: 18.22921 136.0432
## 193: 29.40000 48.0000
## 194: 20.50000 333.0000
## 195: 41.20000 199.0000
## DPT Immunization HepB3 Immunization Measles Immunization Hospital beds
## 1: 66.00000 66.00000 64.00000 0.4322222
## 2: 99.00000 99.00000 95.00000 3.0523077
## 3: 91.00000 91.00000 80.00000 1.8000000
## 4: 57.00000 53.00000 51.00000 0.8000000
## 5: 95.00000 99.00000 93.00000 2.5806250
## ---
## 191: 89.00000 89.00000 95.00000 2.5941667
## 192: 85.68527 81.30779 84.85473 2.9863456
## 193: 73.00000 73.00000 67.00000 0.6672222
## 194: 88.00000 88.00000 93.00000 1.9666667
## 195: 90.00000 90.00000 85.00000 2.3500000
## Basic sanitation services Tuberculosis treatment Urban population
## 1: 49.00617 91.00000 25.754
## 2: 99.18307 88.00000 61.229
## 3: 86.13850 86.00000 73.189
## 4: 51.39362 69.00000 66.177
## 5: 85.52290 72.26667 24.506
## ---
## 191: 87.70814 91.00000 36.628
## 192: 99.37128 76.72564 95.832
## 193: 53.62694 85.00000 37.273
## 194: 31.48123 89.00000 44.072
## 195: 35.77434 84.00000 32.210
## Rural population Non-communicable Mortality Sucide Rate
## 1: 74.246 36.20000 3.60000
## 2: 38.771 6.00000 2.70000
## 3: 26.811 12.80000 1.80000
## 4: 33.823 19.40000 2.30000
## 5: 75.494 17.60000 0.80000
## ---
## 191: 63.372 13.70000 4.70000
## 192: 4.168 22.09954 10.61922
## 193: 62.727 24.80000 4.60000
## 194: 55.928 20.60000 2.70000
## 195: 67.790 27.10000 8.80000
plot_ly(
led,
x = ~ `Tuberculosis Incidence`,
y = ~ `Sucide Rate`,
size = 1,
color = ~ continent,
colors = c("#DF536B", "#61D04F", "#2297E6", "#28E2E5", "#CD0BBC"),
alpha = 0.5,
width = 1.5,
height = 1.5
) %>%
layout(
title = "Распределение частоты самоубийств в зависимости от\nвстречаемости туберкулеза у жителей различных континентов",
xaxis = list(title = "Встречаемость туберкулеза",
zeroline = FALSE),
yaxis = list(title = "Частота самоубийств",
zeroline = FALSE)
)
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Проверим, различаются ли средняя ожидаемая продолжительность жизни в странах Африки и Америки.
Для сравнения применим t-тест Стьюдента. Чтобы понять, необходимо ли нам применять тест Уэлча, оценим равенство дисперсий.
leveneTest(`Life expectancy` ~ continent, data = led %>% filter(continent %in% c("Africa", "Americas")))
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 7.5392 0.007319 **
## 88
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Дисперсии негомогенны. Применим t-тест с поправкой Уэлча (t-тест по умолчанию в R)
for_vis_tt <- t_test(`Life expectancy` ~ continent, data = led %>% filter(continent %in% c("Africa", "Americas"))) %>%
add_xy_position(x = "continent")
И визуализируем его:
led %>% filter(continent %in% c("Africa", "Americas")) %>%
ggplot(aes(x=continent, y = `Life expectancy`))+
geom_boxplot(fill="green")+
labs(
x = "Континент",
y = "Ожидаемая\nпродолжительность жизни",
title = "Распределение ожидаемой продолжительности жизни\nв странах Африки и Америки"
)+
theme(
axis.text.x = element_text(size=16),
axis.text.y = element_text(size=16),
axis.title.x = element_text(size=19),
axis.title.y = element_text(size=19),
plot.title = element_text(size=23, hjust=5),
legend.title = element_text(size=19),
legend.text = element_text(size=16),
legend.position = "bottom",
legend.key.size = unit(1.5, "line"),
strip.text.x = element_text(size=16)
)+
stat_pvalue_manual(for_vis_tt, tip.length = 0)
led_num <- led %>%
select(where(is.numeric) & !Year)
# Шкалируем
led_num_scaled <- scale(led_num)
autoplot(led_num_scaled)
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
Читается этот хитмэп не очень хорошо. Попробуем сделать корплот
# Спирменовская корреляция -- потому что мы не проверяли распределение данных, а наверняка там найдется что-то ненормально-нелинейное
led_num_corr <- corr.test(led_num , method="spearman")
corr_plot <- ggcorrplot(led_num_corr$r,
type = "lower",
outline.col = "white",
lab=TRUE,
method = "circle",
p.mat = led_num_corr$p,
title = "Скоррелированность количественных данных\nв датасете по ожидаемой продолжительности жизни",
legend.title = "Коэффициент\nкорреляции")
corr_plot
# Делаем матрицу дистанций
led_num_dist <- dist(led_num_scaled,
method = "euclidean")
led_num_hc <- hclust(d = led_num_dist,
method = "ward.D2")
fviz_dend(led_num_hc,
cex = 0.1)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
pheatmap(led_num_scaled,
show_rownames = FALSE,
clustering_distance_rows = led_num_dist,
clustering_method = "ward.D2",
cutree_rows = 5,
cutree_cols = length(colnames(led_num_scaled)),
angle_col = 45,
main = "Dendrograms for clustering rows and columns with heatmap")
Мы наблюдаем следующие 2 группы скоррелированных переменных: 1) наличие иммунизации от кори, гепатита B и АКДС 2) экономические метрики: ВВП и ВНД
И 2 группы менее скоррелированных: 1) Безработица, встречаемость туберкулеза, лечение туберкулеза, смертность в ДТП, детская смертность и доля сельского населения. 2) Частота самоубийств, количество больничных коек (страшная вырисовывается корреляция), базовые санитарарные услуги, продолжительность жизни, доля городского населения, подушевой доход.
В данных мы можем выделить 4 субкластера, один из которых чрезвычайно мал (второй сверху на вертикальной дендрограмме).
led_full.pca <- prcomp(led_num_scaled) # Это уже отшкалированные данные, дополнительно не шкалируем
summary(led_full.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion 0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion 0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
## PC15 PC16 PC17 PC18 PC19
## Standard deviation 0.34546 0.26941 0.20224 0.06968 0.0000000000000002735
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.0000000000000000000
## Cumulative Proportion 0.99377 0.99759 0.99974 1.00000 1.0000000000000000000
Первые 6 компонент объясняют 80% вариации. Неидеально, но приемлемо. Посмотрим на соответствующий график.
fviz_eig(led_full.pca, addlabels = T, ylim = c(0, 40))
Посмотрим на переменные
fviz_pca_var(led_full.pca, col.var = "contrib")
Заметим, что все так же выделяется группа иммунизации, и две группы “слабо скоррелированных” переменных с фитмэпа.
fviz_pca_var(led_full.pca,
select.var = list(contrib = 5),
col.var = "contrib")
5 переменных, вносящих наибольший вклад в данные – все те же три переменные о проценте иммунизации, продолжительность жизни и младенческая смертность.
Посмотрим, кто вкладывается в избранные 6 главных компонент
fviz_contrib(led_full.pca, choice = "var", axes = 1, top = 24) # 1
fviz_contrib(led_full.pca, choice = "var", axes = 2, top = 24) # 2
fviz_contrib(led_full.pca, choice = "var", axes = 3, top = 24) # 3
fviz_contrib(led_full.pca, choice = "var", axes = 4, top = 24) # 4
fviz_contrib(led_full.pca, choice = "var", axes = 5, top = 24) # 5
fviz_contrib(led_full.pca, choice = "var", axes = 6, top = 24) # 6
В первую наибольший вклад вносят продолжительность жизни, младенческая смертность, санусловия и соблюдение гигиены и правил приготовления пищи.
Во вторую – переменные вакцинации.
В третью – подавляющий вклад со стороны экономических показателей.
В четвертую – частота самоубийств.
В пятую основной вклад вносит безработица.
В шестую – она же, плюс сопутствующие факторы.
Построим би-2плот
ggbiplot::ggbiplot(led_full.pca,
scale=0,
alpha = 0.5) +
theme_minimal()